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We describe a method for evolving tiie projected Gross-Pitaevskii equation (PGPE) for an interacting Bose gas 
in a harmonic oscillator potential, with the inclusion of a long-range dipolar interaction. The central difficulty 
in solving this equation is the requirement that the field is restricted to a small set of prescribed modes that 
constitute the low energy c-field region of the system. We present a scheme, using a Hermite-polynomial based 
spectral representation, that precisely implements this mode restriction and allows an efficient and accurate 
solution of the dipolar PGPE. We introduce a set of auxiliary oscillator states to perform a Fourier transform 
necessary to evaluate the dipolar interaction in reciprocal space. We extensively characterize the accuracy of 
our approach, and derive Ehrenfest equations for the evolution of the angular momentum. 



PACS numbers: 02.60.Cb,03.75.Hh 

I. INTRODUCTION 

The phenomenal recent progress in experimental efforts to 
produce quantum degenerate dipolar gases [1-5] has brought 
these systems to the forefront of atomic and condensed mat- 
ter physics, driven by a broad range of exciting applications 
[6-13]. Although extensive work has been done on theory 
for the T = dipolar system (e.g. see [10, 14-22]), a gen- 
eral finite temperature theory has yet to be established. The 
long-range character of the dipole-dipole interaction has made 
the development of finite temperature methods more challeng- 
ing. For example, meanfield treatments (which have served as 
the workhorse theory for Bose gases with short-range inter- 
actions) have only been applied to the dipolar gas with ad- 
ditional approximations made to the treatment of exchange 
interactions [23], and quantum Monte Carlo calculations are 
limited to small numbers of particles [24] . 

Recently various classical field methods have become pop- 
ular in the description of ultra-cold Bose gases interacting 
with short range interactions [25-31]. The appeal of these 
methods is that the dynamics of the modes are treated non- 
perturbatively so that non-equilibrium situations or strongly 
fluctuating equilibrium systems (e.g. see [32]) can be accu- 
rately simulated. In Ref. [33] we have developed a quantita- 
tive classical field formalism referred to as c-field theory [34], 
for which the projected Gross Pitaevskii equation (PGPE) is 
the underlying equation of motion. This approach has found 
good agreement with experiment in the critical region of the 
condensation transition [32], and has seen numerous appli- 
cations to regimes where traditional meanfield methods are 
inapplicable (e.g. see [35, 36 1). A key component of c-field 
theory (and the primary distinction from other finite tempera- 
ture classical field theories [27]) that enables it to be applied 
to the quantitative description of experiments is the use of a 
projector, i.e. the explicit restriction of our description to the 



low energy modes of the system. 

In the literature various numerical techniques have been de- 
veloped for for solving the (T = 0) dipolar Gross-Pitaevskii 
equation, such as Crank-Nicholson [21], Fourier pseudospec- 
tral [37], split-operator Fourier transform [14] and split-step 
Fowier transform [22, 38] methods. Underlying all of these 
approaches is the use of a uniform spacial grid which enables 
the efficient evaluation of the dipolar term with Fast Fourier 
transforms. For accurate simulation of 3D dipolar gases these 
approaches require ~ 10^ spatial grid points. In finite tem- 
perature applications the number of grid points corresponds 
to the number of modes that are thermally accessible, and the 
aforementioned approaches tend to have orders of magnitude 
too many modes. Indeed, for typical experimental situations 
of the order of a few thousand modes are appropriate to be 
described by the PGPE [39]. In previous work [40] we have 
found that a practical way to enforce this restriction is by us- 
ing a numerical approach based on a spectral representation 
[41, 42]. 

In this paper we develop the numerical underpinnings of a 
c-field theory for the dipolar Bose gas by introducing a suit- 
able spectral technique for solving the dipolar PGPE. The out- 
line of this paper is as follows. In Sec. n we discuss the dipolar 
PGPE and the spectral representation necessary to implement 
the explicit projection. In Sec. Ill we briefly review the PGPE 
algorithm for the trapped Bose gas with contact interactions, 
before presenting our extension to the dipolar case in Sec. IV. 
In Sec. V we present results characterizing the accuracy of our 
scheme, making comparison to some exactly known matrix el- 
ements and other results in the literature. We also examine the 
convergence of our calculations of equiUbrium properties to 
provide evidence that the scheme we have developed is suit- 
able to making rehable physical predictions. 
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FIG. 1: (a) Schematic diagram sliowing the c-field (C) and inco- 
herent (I) regions of the single particle spectrum for a harmonically 
trapped Bose gas. The energy ecut is usually chosen so that the aver- 
age number of particles in the modes at the cutoff is ricut ~ 1. (b) A 
typical example of an instantaneous c-field density slice for a dipolar 
matterwave with ecut = 23. 



of the c-field technique, and phenomenologically motivates 
the replacement of the quantum field operator for these modes 
by a classical field, i.e. i/jq V'c- This replacement can be 
rigorously justified via a Wigner representation of the many- 
body density matrix, e.g. see Ref. [34]. However, an im- 
mediate consequence of this development is that the c-field 
formalism must be restricted to the low energy modes of the 
system where this field replacement is valid (i.e. the c-field re- 
gion, C, shown schematically in Fig. 1(a)). To formalize this 
restriction we introduce a projector, Vc 



rc{F{^)} ^ }_^M^) / d'xV:(x')^^(x'), (3) 

C ^ {n : e„ < Ecut}, (4) 
where (j)n{^) and e„ are eigenstates of Hq, i.e. 

e„0„(x) = _ffo0n(x), (5) 

and the (single particle) energy cutoff, e^ut. is the single pa- 
rameter we use to define the c-field region [43]. The action 
of Vc in Eq. (3) is thus to project the arbitrary function F(x) 
into the c-field region. 

The equation of motion for the c-field treatment of a Bose 
gas is the projected Gross-Pitaevskii equation (PGPE). For the 
case of a gas of particles interacting via short range and long 
range dipole interactions, the PGPE takes the dimensionless 
form 



+ / dVl/i5(x-x')|^Ac(x')rV'c(x) k (6) 



where 



D 



1 - 3 cos2 



(7) 



II. FORMALISM: DIPOLAR PGPE 

Our interest is in a system of bosonic particles confined in 
a harmonic potential, described by the single particle Hamil- 
tonian 



Ho 



V^o(x) 



1 



V2 + T/o(x), 



1 ^ 



2 2 



(1) 



(2) 



where Xj — LOjjuj is the relative trap frequency in each di- 
rection i = {x, y, z}. To obtain this dimensionless form we 
have used harmonic oscillator units of length xq — ^hjmuj, 
energy Eq = huj and time to = l/uj, with m the particle mass 
and uj a convenient reference frequency. 

Near thermodynamic equilibrium the low energy modes of 
the system are highly occupied and their dynamics are domi- 
nated by classical fluctuations. This observation is at the heart 



is the dipole interaction potential with r = x'^ + + 
and 9 the angle between x and the z axis (the axis along 
which the dipoles are polarized). Here we have introduced 
the dimensionless s-wave (contact) interaction parameter C — 
AiraNQ/xQ, with a the s-wave scattering length, and the di- 
mensionless dipole interaction parameter D = NccPrn/h^xo, 
with d the dipole moment. For convenience we take the field 
ijjc to be normalized to unity so that the number of c-field 
atoms, Nq, appears explicitly in the interaction parameters. 

The usual strategy for dealing with the dipolar interaction 
is to make use of the Fourier transformed density and dipolar 
interaction potential 



(8) 
(9) 
(10) 



n(k) = J d-'xe-'^-^\M^)\', 
^^(k) = J d'xe-^'^-^VDi^), 

= -l^[l-3cos20fe], 
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where 6k is the angle between k and the kz axis. Thus making 
use of the convolution theorem we have 

$(X) = y(iVl/z3(x-x')|Vc(x')|', (11) 



d-'fce'''-^VD(k)n(k). 



(12) 



The main concern of this paper is to develop a suitable method 
for evaluating $(x) in a manner appropriate for use in the 
PGPE formalism. We emphasize that the modes of the sys- 
tem are of central importance in the PGPE and care must be 
taken in numerical implementations to ensure the modes are 
faithfully represented. This point is made clear with reference 
to Fig. 1(b), which shows a snapshot of the c-field density and 
reveals the appreciable occupation of every mode in the c-field 
region. 

We also note that the energy functional for the dipolar 
PGPE is 

E[i,c] = J d'xrcHoi^c+\ J d^xC\M^)\^ (13) 



which forms an important constant of motion for the system. 
In a similar manner to how we dealt with the dipolar part of 
the PGPE, it is convenient to evaluate the dipolar energy term 
in Fourier space as 



d3a;$(x)|^c(x)|' 



^11 



d3fcVD(k)n(k)fi(-k). 



(14) 



A. Spectral representation 



It is most convenient to expand the c-field in a spectral basis 
of the single particle states, i.e. 



V'c(x,i) = ^Cn{t) ?!>„(x), 



(15) 



nec 



where the {c„} are complex amplitudes. The projection is 
explicitly implemented by limiting the summation indices in 
(15) to the set of values specified in Eq. (4) defining the c-field 
region. 

B. Mode evolution 

Having used the modes of Hq as the spectral basis and to re- 
aUze the projector, we follow the Galerkin approach (i.e. pro- 
jecting Eq. (6) on to our spectral basis) to obtain the evolution 
equation for the mode amphtudes 



dCn 

dt 



-i [CnCn + Gn] 



(16) 



where 



Gn = J rf'x <^;(X) [C\i,cf + i>cM, (17) 



is the nonlinear matrix element. Once these nonlinear ma- 
trix elements are evaluated, the evolution of the system can be 
calculated using numerical algorithms for systems of ordinary 
differential equations, e.g. the Runge-Kutta algorithm. Since 
this is a well-understood area of numerical mathematics we 
do not concern ourselves with the details of the propagation 
algorithm, but instead focus on evaluating Eq. (17). 

In principle the nonhnear matrix elements between spectral 
basis functions can be computed exactly. Defining 



^npqr 



d3a:C(x) [C0;(x)0,(x),/),(x)+ (18) 
J dPx' FD(x-x')(/);(x')</)g(x')<^.(x)] , 

which can be calculated analytically (c.f Appendix C), and 
expanding the c-field in terms of its spectral representation 
we see that 



Gn — ^ ^ Inpqr^p^q^v 

{p,q,r}ec 



(19) 



While being exact, evaluating this expression is prohibitively 
slow, requiring 0{M'^) operations, where M is the number 
of modes in the c-field region. In contrast, the approach we 
develop here is 0(M^/^), and thus suitable for simulating real 
systems in a reasonable amount of time (e.g. simulations of 
the order of hours to days on a commodity PC). 



C. Separability 

In what follows we take the trap to be isotropic, and set all 
\j = 1, for simphcity of notation [44]. An important fea- 
ture of the basis states (i.e. eigenstates of Hq) is that they are 
separable into ID eigenstates, i.e. 



<^n(x) ^ lfa{x)^p{y)'p^{z), 



(20) 
(21) 
(22) 



where {ipa{x)} are eigenstates of the ID harmonic oscillator 
Hamiltonian, i.e. 



I <f 1 2 



iPa{x) = ea(Pa{x), 



(23) 



with eigenvalue = + |), for a a non-negative integer. 

For clarity we use greek subscripts to label the ID eigen- 
states, so that the specification of the c-field region in (4) be- 
comes 



C = {a, f3,j : Sa + + < Ccut}- 



(24) 



Within the c-field region there exists (« Ccut) distinct ID 
eigenstates (i.e. (pa) in each direction, and thus 
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(25) 



3D basis states (</>„) in the c-field region. 
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III. REVIEW OF STANDARD PGPE ALGORITHM 

We first begin by reviewing the PGPE algorithm we have 
developed for the case of local interactions. This algorithm 
uses Gauss-Hermite quadrature to evaluate the (local) nonlin- 
ear term exactly in an efficient manner. For a complete ac- 
count we refer the reader to Ref. [40]. 

A. Evaluating the matrix elements 

To begin we note the harmonic oscillator states are of the 
form 

iPo^ix) = haH„{x)e-''"/^, (26) 

where = [2"a!v^~^/^ is a normaUzation constant, and 
Ha (x) is a Hermite polynomial of degree a, defined by the 
recurrence relation 

Hc,+i{x) = 2xHa{x) - 2aHa-i{x), a = 1,2,..., (27) 

with Ha{x) = 1 and Hi{x) = 2x. 

Thus, the field (at any instant of time) can be written as 

Vc(x, t) = Q{x, y, z)e-^-"+y"+'"^/^ (28) 

where 

Q{x,y,z)= ^ Cai3^{t)haHa{x)hf3Hf3{y)h^Hj{z), 

{"/97}ec 

(29) 

is a polynomial that, as a result of the cutoff, is of maximum 
degree — 1 in the independent variables. 

Similarly, it follows that because the interaction term (17) 
is fourth order in the field, it can be written in the form 

Ga0j = J d^x e-2(^'+^'+^')Pa^^(x, y, z), (30) 

where 

Pai3j{x,y,z) = C haHc{x)h0Hi}{y)h^H^{z) 

x\Q{x,y,z)\'^Q{x,y,z), (31) 

is a polynomial of maximum degree 4 (M^ — 1) in the in- 
dependent variables. To evaluate these integrals, we note the 
general form of the Nq point Gauss-Hermite quadrature 

/ dxw{x)f{x) »^Wjf{xj), (32) 

where is a Gaussian weight function, and the Nq values 
of Wj and xj are the quadrature weights and roots, respec- 
tively. This quadrature is exact if f{x) is a polynomial of 
maximum degree 2Nq — 1. 

Identifying the exponential term in (30) as the weight func- 
tion for quadrature, the integral can be exactly evaluated us- 
ing a three-dimensional spatial grid of 8 {M^ — 1)^ points (i.e. 
2 {Mx — 1) points in each direction [45]), i.e. 

Gafij = y^,WiWjWkPa0j{Xi,Xj,Xk), (33) 
ijk 



where Xi and Wi are the 2 {M^ — 1) roots and weights of the 
ID Gauss-Hermite quadrature with weight function w{x) = 
exp(— 2x^). Note, that the isotropy of the trapping potential 
(for the numerical examples considered in this paper) results 
in identical quadrature grids in all spatial directions in our ex- 
ample. 



B. Overview of the numerical algorithm 

Here we briefiy overview how the quadrature described 
above can be efficiently implemented numerically. We require 
the transformation matrices, given by ID basis states evalu- 
ated on the quadrature grid, i.e. 

Uia = V>a{Xi), (34) 

to be pre-calculated. Because the transformations are block 
diagonal, i.e. applied across the directions independently at 
computational cost O(M^) = ©(M'*/^) (see Eq. (25)), we 
will make use of the simphfying notation 

^ UtaUji3Uk^Caf3^{t) ^^^Usa-Ca-, (35) 
{a/37} e C cr 

where cr = {aPj} and s = {ijk}, and it is understood that 
c<T = Ca/37, and C/s<T = Uial/jfjUk-y. 

Starting from the basis set representation of the field (i.e. 
{cap-y}) at an instant of time t, the steps for calculating the 
matrix elements are as follows: 

Step 1: Transform from spectral to spatial representation: 

V'c(Xs) = ^J/so-Co-, (36) 

(T 

where Xg = {xi,Xj,Xk). 

Step 2: The quadrature integrand of the nonlinear matrix ele- 
ment (17) is constructed by appropriately dividing by 
the weight function and pre-multiplying by the weights 
[46], i.e. 

5(x«) = Wse''\^'\"c\M^s,t)\''M^s,t), (37) 

where Ws = WiWjWk- 

Step 3: The inverse transform of g{^s) yields the desired ma- 
trix elements: 

= Y.^t^9{^s)- (38) 



The slowest step in this procedure is carrying out the basis 
transformation (steps 1 and 3), which requires 0{M^/^) float- 
ing point operations when carried out as a series of matrix 
multiplications. Thus, the overall algorithm is 0{M^/^). 
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IV. EXTENSION TO CALCULATE THE DIPOLAR TERM 

To treat the dipolar term we need to augment step 2 in the 
standard harmonic PGPE algorithm (see Sec. IIIB). To do 
this we want to Fourier transform the density associated with 
V'c(xs, to form Eq. (12). It is not convenient to use a fast 
Fourier transform because ?/'c(xs, t) is evaluated on a nonuni- 
form grid (i.e. quadrature grid). Interpolation to a uniform 
grid would be computationally expensive and would introduce 
a source of considerable error, especially since the quadrature 
grids tend to be quite sparse (see discussion in Sec. IV B). 

Here we show how an auxiliary harmonic oscillator basis 
can be used to perform the Fourier transform exactly. Follow- 
ing similar arguments to those made in Sec. Ill A, the c-field 
density, n(x) = |V'c(x)|^, is of the form 



n(x) = R{x,y,z)e 



(39) 



where i? is a polynomial of maximum degree 2 (Mr — 1) in 
the independent variables. 

Introducing a set of auxiliary harmonic oscillator states. 



(40) 



which differ from the spectral basis oscillator states by a factor 
of 2 in the argument of the exponential (chosen to match the 
exponential part of Eq. (39)). Indeed, these states are eigen- 
states of the operator 



+ 2x' 



(41) 



i.e. harmonic oscillator with twice-as-tight trapping potential, 

and expressions for and Ha (x) can be obtained by noting 
that these modes relate to the usual dimensionless oscillators 
by a simple scaling Xaix) = 2^/'^ipce{\/2x). 

The auxiUary oscillator states form an orthonormal basis, 
and because of their appropriate exponential factor, we can 
exactly represent the density (39) as 



n(x) 



ic?<TXcr(x) 



(42) 



where d^r <-> dajs^ is a set of 8M^ real coefficients, with 
Xo-(x) = Xa{x)x0{y)Xi{z)- Indeed, because the {x<t} are 
an orthonormal basis, we have that 



d. 



= j d3a;e-2|-l's<,(x), 



(43) 
(44) 



where in the second line we have collected exponential and 
polynomial terms separately, with 



5<,(x) = e2|-l\;(x)n(x), 



(45) 



a polynomial of degree 4(Mj. — 1) in the independent vari- 
ables. Thus the integration (44), like that in Eq. (30), has same 



weight function and maximum degree of polynomial order. 
Thus Eq. (44) can be calculated exactly with the same quadra- 
ture (i.e. roots {xi} and weights {wj}) as used in Eq. (33), i.e. 



da 



Ws S„{Xs). 



(46) 



The harmonic oscillator states are eigenstates of the Fourier 
transform operator with eigenvalue (— i)", i.e. 



1 



(27r)i/2 



dxe 



Xoc{x). 



(47) 



Thus knowledge of the basis amplitudes dap-^ allows us to 
efficiently and precisely construct the Fourier transform of the 
classical field density, i.e. 



n(k) = Y.d,{-i)-\\'^H.{n 



(48) 



where \cr\i — a + /3 + 7 is the one norm of cr, noting 
that {a, /3, 7} are non-negative. We can now construct the 
integrand of the dipolar interaction term in Fourier space, 
i.e. yD(k)n(k) appearing in Eq. (12), which needs to be in- 
verse Fourier transformed to obtain ^(x). This can be done 
using the inverse of the procedure we used to obtain n(k), i.e. 
via the expansion of $(x) in the auxiUary oscillator states 



(49) 



where 



f^ = yd3/e(i)-Hli^;(k)Vb(k)n(k). (50) 

Expression (49) is approximate because Vd{^) is not of the 
form of a finite-degree polynomial, and thus yD(k)n(k) can- 
not be represented exactly in the oscillator basis - an approx- 
imation we investigate in Sec. V. 

To numerically evaluate the fap^y we again make use of a 
Hermite-Gauss quadrature with roots {A;,} and weights {tDj}, 



I.e. 



where 



= Y^^MK), (51) 



T^(k) = e2|k|'x;(k)T>D(k)n(k). (52) 



Note the number of A;-grid quadrature points is in principle 
arbitrary, but should be at least 2Mx in each direction. We 
can use the number of points to control the accuracy of the 
matrix element. 



A. Spectral dipolar algorithm summary 

Step 1: Transform from spectral to spatial representation: 

V'c(Xs) = ^J/so-Co-. (53) 
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Step 2a: The weighted position density is constructed 

/(x,) = u;3e^l-»l'|^c(x«)|'. (54) 
Step 2b: We compute the Fourier transformed density as 

n(kt) = Y^W^tfi^s), (55) 

where t = {uvw} are the indices which label the 
Fourier space grid points. Here we have introduced the 
pre-computed transformation matrix, 

Wir = $^(-i)"Xa(fc.)Xa(a;i). (56) 

a 

which combines both steps of the Fourier transform into 
one (i.e. n(x) — > da^j and da/3-y Ti(k)). 

Step 2c: The product with the dipole interaction potential is then 
formed in Fourier space 

/(kt) = %e2|''*l'y^(kt)n(kt). (57) 

[Or with the replacement VoCkt) ^oi^t)' ^ 
rected dipolar interaction, as discussed in Sec. VB]. 

Step 2d: Inverse transforming yields 

$(X3) = ^W:J{kt). (58) 
t 

Step 2e: Short range and dipolar interaction terms are then com- 
bined into a single integrand 

ff(x,) = «;,e2|-»l' [C|Vc(x3)|' + $(x«)] Vc(x.).(59) 

Step 3: Inverse transforming this integrand yields the desired 
matrix elements: 

= ^C/;^5(X3). (60) 

s 

Steps 1, 2b, 2d, and 3 are 0(M^/^). Since the algorithm 
involves twice as many transformations as the non-dipolar 
PGPE case, each evaluation of the (and hence each time 
step) takes approximately twice as long. 

B. Possibility of using fast Fourier transformations 

Having presented our spectral algorithm we are now able 
to comment on the alternative procedure of computing $ us- 
ing fast Fourier transformations (FFTs). To do this requires 
several modifications to the algorithm, which we briefly sum- 
marize. In step one, in addition to computing V^c(xs) on the 
quadrature grid for the short range interaction, we will need 
a new transformation Us to obtain tpci^s) on the uniformly 
spaced grid {Xg}. Following standard procedures (e.g. see 
[14]) we can then obtain $(xs) using two FFTs. This step 



is more efficient than our procedure using Wst in the spectral 
algorithm, but we will likely require more x^ grid points for 
the Fourier representation to provide an adequate representa- 
tion of trapped field. Additionally, the efficiency of the FFTs 
is offset by the need to interpolate Xs back onto a quadra- 
ture grid for step 3 (if performed on a uniform grid this last 
step is highly inaccurate without a prohibitively large num- 
ber of points). Due to the added complexity of the FFT algo- 
rithm, and that approximations occur in the algorithm at sev- 
eral places, we have decided not to investigate this any further 
in this work. 



V. ACCURACY OF APPROACH 

Step 2d of our numerical algorithm for the dipolar PGPE 
is approximate and requires investigation to justify that it 
is sufficiently accurate to be useful. The PGPE formalism 
places strong constraints on the underlying numerical algo- 
rithm which restrict how we might improve the accuracy. In 
particular, the c-field region is defined by Ccut. and hence 
is dictated by the physical system under consideration (i.e. 
temperature, number of atoms) and is not a parameter that 
can be arbitrarily varied. Instead, for fixed M^, we would hke 
to understand: (i) The accuracy of the matrix elements 
(ii) What ways we have for controlling this accuracy? (iii) 
What level of accuracy is needed for making reliable physical 
predictions? 

Here we investigate two methods of improving the accuracy 
of the matrix elements. The first method, which we discuss in 
Sec. V A, is by increasing the order of the fc-space quadrature. 
The second method is to use a modified (finite range) interac- 
tion potential, which we present in Sec. VB. We then charac- 
terize the effect of these adjustments using various tests. We 
finally tum to addressing what level of accuracy is required to 
make useful predictions with the PGPE theory. 

A. Fourier quadrature grid 

The two quadrature grids {xg} and {kt} are central to the 
computation of the nonlinear matrix elements in our algo- 
rithm. Since the weight functions are known for each quadra- 
ture they are completely specified by the number of points, 
i.e. the parameters 

N^. The number of quadrature points along each direction 
in the position space x grid. 

Nk'. The number of quadrature points along each direction 
in the Fourier space k grid. 

First, we note that for given ecut (i e. M^) the transform to 

fc-space is exactly invertible (i.e. n(x) n(k) n{x)) if 
we choose > N^, Nk > N^, where we have defined the 
reference values 

= 2M^-1, (61) 
N° = (62) 
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Note, that the invertible requirement is met with = 2Mx — 
1, but we choose 2Mx to avoid having an odd number of 
points which ensures that there is no quadrature point at k = 
where Vd is singular. 

With the inclusion of the dipolar potential it is beneficial to 
increase the number of momentum grid points beyond A^;^ to 
obtain better accuracy for step 2d. In the results we present 
below we will indicate the increase in momentum grid points 
over the reference value as AA^^, i.e. 

Nk = N^ + ANk. (63) 

We do not alter from the reference value N^, as this has 
no effect on the accuracy of the algorithm. 

B. Corrected dipolar interaction 

Ronen et al. [38] have demonstrated a useful procedure for 
improving the convergence of the numerical evaluation of the 
dipolar term for low energy states in Bogoliubov calculations. 
They noted that the poor convergence of this term arises be- 
cause the Fourier transformed interaction, VoCk), is singular 
at the origin (where r;,(k) is typically large) due to the long 
range character of the interaction. Ronen et al. suggested 
the use of the Fourier transform of the dipolar interaction re- 
stricted to a spherical domain of size R, i.e. the Fourier trans- 
form of 

f D{l-Scos^9)/r^ r<R, 
y^(x) = I (64) 
0, otherwise. 

This has the analytic transform 

V3i^) = f (l + 3^ -3^) (3cos^..-l), 

(65) 

which we shall refer to as the corrected dipolar interaction, 
having the feature that it is less rapidly varying near k = 0. 
This approach seems reasonable as we are studying a trapped 
system of finite spatial extent, and thus the sharp behavior of 
the uncorrected potential (Vd) at k = 0, arising from inter- 
actions over long length scales, cannot be physically relevant. 
Ronen et al. justify using V§ as it prevents the "long range 
interactions between copies of condensates" arising from the 
periodicity of their Fourier based calculations. 

More generally, the use of V§; can be justified by noting 
that sharp features in the interaction potential are not accu- 
rately calculated on a finite quadrature grid (or Fourier grid). 
In practice if these sharp features are left in the numerical cal- 
culations they are misrepresented by the finite quadrature and 
interfere with lower order matrix elements (often referred to 
as aliasing in the Fourier case), leading to their slow conver- 
gence as the number of quadrature points is increased. 

Choice for R 

An immediate issue to investigate is the optimal choice of 
the length scale R. For the our trapped system the character- 



istic size is given by the classical turning point Itp ~ \/2M^ 
(in computational units), since ecut ~ Mr- 

To investigate the accuracy of our algorithm as we vary R 
used in the corrected dipolar interaction we consider the pure 
dipolar matrix element (for D = 1): 

Zl= j d?xSx' (j)*^{^)VD{-i^--s.')\(t>u{-^)?M-^)- (66) 

to be distinguished from the general matrix element which re- 
quires four distinct oscillator state labels. In practice we eval- 
uate this as follows: we take — 6^^^, and then compute 
the nonlinear matrix elements Gt using our algorithm (see 
Sec. IV A), and identify = G-r- These pure dipole ma- 
trix elements are useful for characterizing the accuracy of the 
algorithm and we will make use of these in several apphca- 
tions. It is convenient to indicate the matrix elements under 
consideration using the notation Z^^^, as established earlier 
(E.g. see Eq. (20)). 

For the purposes of studying the dependence on R we will 
consider four non-trivial matrix elements shown in Table I for 
which we calculate the values exactly using an analytic ap- 
proach discussed in Appendix C. The results of our algorithm 
are shown in Fig. 2 and confirm that the corrected interaction 
potential, V^, has considerable advantage over the bare po- 
tential, Vd, for certain values of R. However, it is clear that 
there is quite complex structure in the variation of with R 
and there is no single value of R for which all matrix elements 
obtain the smallest relative error. 



^000 = 0.037612 ••• 

^000 = 0.009498 ••• 

ZiSo -0.003 908 ■■• 

Ztti 0.000 462 291 •■• 



TABLE I: Values of pure dipole matrix elements considered (see 
text) 

To interpret these results it is useful to qualitatively classify 
the matrix elements into two categories: 
Low order matrix elements: These are matrix elements that 
involve low order oscillator states, i.e. those with quantum 
numbers much less than Af,,. (i.e. the cases in Figs. 2(a) and 
(b)). For these cases the typical density variations are well- 
resolved on the quadrature grids and for both cases we see 
that R w ^/2Mx is the optimal value for obtaining such ma- 
trix elements with small relative error. This value appears to 
be universally good for low order matrix elements 
High order matrix elements: These are matrix elements that 
involve oscillator states with quantum numbers comparable 
to Mr (i.e. the case in Figs. 2(c) and (d)). For these cases 
the typical density variations are rapid on the quadrature grids 
and R « a/2M^ is clearly not the optimal value for obtaining 
such matrix elements with small relative error. The location of 
the minimum relative error (e.g. R « O.A\/2Mx in Fig. 2(c), 
R « 0.65\/2MjJ^ in Fig. 2(d)) appears to vary appreciably 
with the particular high order matrix element, so that there is 
no universally good value. 
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FIG. 2: Relative error in the pure dipole matrix elements (a) Zqoq, 
(b) Z^SS, (c) Ziti and (d) as 7? is varied. Results: (solid line) 
calculated using VoiM) ; (circles) calculated using V^{k). Other 
parameters are at the reference values with = 16. 



low order matrix elements is i? w ^/2Nk- This can be un- 
derstood as follows: the use of the corrected interaction in- 
troduces an infrared cutoff in Fourier space at the wavevector 
scale fccut ~ 1 /R, which for the reference case Nh = is 
approximately equal to the spacing between k grid points near 
A; = 0. As we increase the k grid resolution improves, 
i.e. smaller wavevectors are resolved and a longer R value is 
needed to represent the correspond longer wavelengths. 
C. Energy convergence for a Gaussian density 

Ronen et al. [38J have checked the accuracy of their numer- 
ics by evaluating the dipolar energy functional [47] 

Id = j j (fx(fx'VD{x-x')n{x')n{x), (67) 

for the case of D = 1 and the Gaussian density 

^i-) = ^2^-'''^^'''''-'^ ^ (68) 

for which the exact result is 

Jz) = 0.038 670 861 ••• . (69) 

The results of Ronen et al. are shown in Table 11, and clearly 
reveal the large improvement they obtained by using the cor- 
rected dipolar interaction. 

It is not possible to directly compare our harmonic oscilla- 
tor approach since we do not have independent control of the 
spatial extent and number of grid points. However, we can 
vary and check convergence [48]. For this case we use 
an isotropic harmonic oscillator potential, so that the density 
(68) cannot be simply related to any finite superposition of 
eigenmodes of Hq. Thus, we expUcitly construct n(x) on the 
quadrature grid before performing the normal transformations 
to make ^'(xs). To calculate the energy functional we then 
evaluate 



In what follows we will take R = -y/2M^. We make this 
choice because this appears to universally improve the accu- 
racy of the low order matrix elements by at least several orders 
of magnitude over the imcorrected values, while only having 
a minor detrimental effect on the accuracy of the higher order 
matrix elements. The cases presented in this section have been 
for the reference value (N^ = 2Mx) of Fourier grid points. If 
additional Fourier points are added the best R value for the 



/D = 5]^se'l"^l'$(xs)n(xs). (70) 

s 

The results shown in Table 11 reveal quaUtatively similar be- 
havior to those observed in by Ronen et al. , i.e. we see that the 
accuracy of the calculation improves gradually as the number 
of points increases, and a rather dramatic improvement in the 
accuracy if the corrected dipolar interaction is used. 



D. Pure dipole matrix element convergence 

In this section we investigate the effect of increasing the 
number of k grid points on the accuracy of pure dipolar ma- 
trix elements. Typical results for the relative error are shown 
in Fig. 3, with the corresponding exact matrix element val- 
ues given in Table 1. Figures 3(a) and (b) show the charac- 



teristic behavior for the low order matrix elements, indicat- 
ing the general trend that these matrix elements improve con- 
siderably with AiVfe. For the higher order matrix elements 
[see Figs. 3(c) and (d)], the improvement in the relative error 
is much more gradual, but quite significant considering the 
rather low relative accuracy of these matrix elements in the 
reference configuration. The case seen in Fig. 3(d) shows that 
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Results of Ref. [38]: 


7? = 8,iV = 32 


Relative error 

R^8,N = 64 R=16,N = 64 


i?=16,7V = 128 


Using V'z,(k) 
Using t/«(k) 


2.7 X 10"^ 
-1.1 X 10"^ 


2.7 X 10"^ 
-1.1 X 10"^ 


8.6 X 10"^ 
1.8 X 10"* 


8.6 X lO"'^ 
-4.4 X 10"" 


Our results: 


= 16 


Relative error 
=32 


=64 




Using Vi5(k) 
Using V«(k) 


-1.7 X 10"^ 
-2.9 X 10"^ 


-3.1 X 10"^ 
-1.9 X 10"^ 


-5.5 X 10"* 
8.3 X 10"^ 





TABLE II: Relative error of the dipole interaction energy. We compare to the results of Ronen et al. [38] using a 3D FFT method on a cubic 
grid of extent [—R, R\ with AT points in each direction 
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by increasing AN^ we can make the error in the corrected 
interaction matrix element smaller than the uncorrected value. 



E. Random state convergence 




FIG. 3: Relative error in the pure dipole matrix elements (a) .^ooO' 
(b) Z'SiM, (c) Zllland (d) Zf^a, as AAffc is varied. Results: (grey 
squares) calculated using Vd{}^) ; (black circles) calculated using 
V§{k). Other parameters: Mj; = 16 and E = \/2Nk (see text). 



FIG. 4: Relative error in the random state matrix elements. Results: 

(grey square) calculated using Vi)(k); (black circles) calculated us- 
ing Voiy^)- Results for a randomized state with (a) Mx = 10 and 
(b) Mx = 30. 

The pure dipole matrix elements considered so far are use- 
ful for tmderstanding the general effects of using ^^(k) and 
changing AiVfe. However, for the purposes of understanding 
the PGPE in operation, a more appropriate test is to determine 
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the nonlinear matrix elements, Go-, for a randomized state Cn-. 
We can then determine the combined effect of altering V^(k) 
and AA^i; by examining 

where Go- refers to the approximate matrix elements, 
refers to the more accurately calculated matrix elements (see 
below), and ||Ao. p = Ylia- \^<^\'^- ^he matrix elements Go- 
determine the transitions between the bare spectral states in 
the PGPE (since iJp is diagonal in that basis) and thus 5G 
measures the extent to which our approximate evaluation of 
Grr matches the more accurate value G^. We note that this 
differs from the earlier consideration of pure matrix elements 
because a large relative error in a small matrix element (typi- 
cally the case for high order modes) has Uttle effect on 5G. 



In Fig. 4 we show results for 5G for cases where G^ is 
evaluated using the bare and corrected dipole interaction, and 
for various values of AA^^.. Our pseudo-random state is repro- 
ducible, with procedure outlined in Appendix A. The accurate 
values, G^, are calculated using our algorithm with A^^ = 128 
quadrature points and . 

The results in Figs. 4(a) and (b) are for = 10 and 
Mx = 30, respectively. In both cases the corrected dipole 
matrix element is more accurate, and converges more rapidly 
with ANf.. For larger the convergence rate is less rapid, 
due to the increase in higher order matrix elements which our 
previous results show to converge more slowly. 



F. Propagation convergence 



Relative Tolerance 




Number of steps 


5N 




SE 




5L, 




sx 




SX' 




10-* 





362 


-2.8x10" 


-4 


2.4x10" 


-3 


4.8x10" 


-2 


1.4x10" 


-3 


8.3x10" 


-2 




10 


346 


-2.9x10" 


-4 


2.5x10" 


-3 


3.2x10" 


-2 


1.6x10" 


-3 


2.2x10" 


-2 




20 


342 


-3.0x10" 


-4 


2.6x10" 


-3 


1.9x10" 


-2 


1.7x10" 


-3 


7.0x10" 


-3 




30 


345 


-2.9x10" 


-4 


2.5x10" 


-3 


1.0x10" 


-2 


1.6x10" 


-3 


2.8x10" 


-3 




40 


348 


-2.9x10" 


-4 


2.5x10" 


-3 


5.2x10" 


-3 


1.6x10" 


-3 


1.8x10" 


-3 


10"^ 





570 


-2.8x10" 


-5 


2.5x10" 


-4 


4.8x10" 


-2 


1.5x10" 


-6 


8.0x10" 


-2 




10 


551 


-2.9x10" 




2.5x10" 


-4 


3.2x10" 


-2 


1.5x10" 


-5 


2.1x10" 


-2 




20 


554 


-2.9x10" 


-5 


2.5x10" 


-4 


1.8x10" 


-2 


1.5x10" 


-5 


5.2x10" 


-3 




30 


554 


-2.9x10" 


-5 


2.5x10" 


-4 


1.0x10" 


-2 


1.5x10" 


-6 


1.1x10" 


-3 




40 


557 


-2.9x10" 


-5 


2.5x10" 


-4 


5.0x10" 


-3 


1.5x10" 


-6 


1.6x10" 


-4 


10"^ 





857 


-2.9x10" 


-6 


2.5x10" 


-5 


4.8x10" 


-2 


1.5x10" 


-7 


8.0x10" 


-2 




10 


868 


-2.9x10" 


-6 


2.5x10" 


-5 


3.2x10" 


-2 


1.5x10" 


-7 


2.0x10" 


-2 




20 


866 


-2.9x10" 


-6 


2.5x10" 


-5 


1.8x10" 


-2 


1.5x10" 


-7 


5.1x10" 


-3 




30 


868 


-2.9x10" 


-6 


2.5x10" 


-5 


1.0x10" 


-2 


1.5x10" 


-7 


1.1x10" 


-3 




40 


872 


-2.9x10" 


-6 


2.5x10" 


-5 


5.0x10" 


-3 


1.5x10" 


-7 


1.4x10" 


-4 


10-^ 





1336 


-2.9x10" 


-7 


2.5x10" 


-6 


4.8x10" 


-2 


1.5x10" 


-9 


8.0x10" 


-2 




10 


1328 


-2.9x10" 


-7 


2.6x10" 


-6 


3.2x10" 


-2 


1.6x10" 


-9 


2.0x10" 


-2 




20 


1344 


-2.9x10" 


-7 


2.6x10" 


-6 


1.8x10" 


-2 


1.6x10" 


-9 


5.1x10" 


-3 




30 


1332 


-2.9x10" 


-7 


2.6x10" 


-6 


1.0x10" 


-2 


1.6x10" 


-9 


i.ixiO" 


-3 




40 


1341 


-2.9 xiO 




2.5 XiO 


(i 


5.0xi() 




1.6 XiO 




1.4 XiO" 


-4 


10"® 





2079 


-2.9x10" 


-8 


2.5x10" 


-7 


4.8x10" 


-2 


1.5x10" 


11 


8.0x10" 


-2 




10 


2089 


-2.9x10" 


-8 


2.6x10" 


-7 


3.2x10" 


-2 


1.5x10" 


11 


2.0x10" 


-2 




20 


2085 


-2.9x10" 


-8 


2.6x10" 


-7 


1.8x10" 


-2 


1.5x10" 


11 


5.1x10" 


-3 




30 


2088 


-2.9x10" 


-8 


2.6x10" 


-7 


1.0x10" 


-2 


1.5x10" 


11 


1.1x10" 


-3 




40 


2092 


-2.9x10" 


-8 


2.6x10" 


-7 


5.0x10" 


-3 


1.5x10" 


11 


1.4x10" 


-4 



TABLE III: Convergence properties of evolution algorithm. The relative error tolerance of the adaptive step Runge-Kutta algorithm, number 
of steps needed to obtain that error tolerance, and the quantitative measures 6N, SE, SL^ and SX are shown (see text). Other parameters: 
T = 1, C = 500, D = 500, ecut = 23 and the initial state is a thermalized state with energy E = 10.0 (see text). All results computed using 
the corrected dipole interaction. 

I 



Here we present some evolution convergence results for our algorithm. We have used an adaptive step Runge-Kutta- 
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Fehlberg algorithm to evolve the dipolar PGPE with a speci- 
fied relative error tolerance. For all the results presented in the 
remainder of this paper we use the corrected dipole interaction 
so as to benefit from its generally more accurate evaluation of 
the matrix elements. Since computing the matrix elements for 
our harmonically trapped algorithm is of computational cost 
0{M^^^) the development of higher order or more efficient 
propagation algorithms would be desirable (e.g. see Refs. [49- 
52]), although we do not address this issue further here. 

We test our algorithm by propagating an initial state for- 
ward in time by an amount T = 1. The system we consider 
has interaction parameters C = 500 and D = 500, and is in 
an isotropic trap potential with energy cut off Ccut = 23, for 
which M = 2024 modes lie in the c-field region. To provide 
a useful analysis of the regime that the PGPE approach is nor- 
mally used, we take an initial state of energy E = 10.0 (as 
given by Eq. (13)) after it has been propagated to thermalize 
for 25 ttap periods. This state has the desirable feature that 
all the modes of the field are appreciably occupied, and thus 
provides a more stringent test of the evolution. 

In Table III we examine the evolution convergence as we 
vary both the integration tolerance and ANk, using the fol- 
lowing measures: 



M 



5N 
SE 

SX 



i-^|c,(r)|2, 

E[V>c(x,r)]-g[V;c(x,0)] 

i?[^c(x,0)] 
(L,(r))-(L,(0)) 



M 



J2\cj{T)-cf{T)\' 



(72) 

(73) 
(74) 

(75) 



i.e. the change in normalization (6N), the relative change in 
energy (SE), the relative change in the z component of an- 
gular momentum (SL^), and a difference measure of the final 
states (6X), where cf{T) are the mode ampUtudes at time T 
of a more accurate simulation (discussed below). The quan- 
tity SX provides a direct test of the field convergence at the 
final time. However, the other quantities considered relate to 
constants of motion, which are useful in practice as they pro- 
vide a characterization of the accuracy without the need for 
running additional simulations. 

a. Normalization The dipolar PGPE formally preserves 
the normalization of the field. Our results in Table 111 show 
that this quantity, as defined in Eq. (72), is dependent on the 
tolerance of evolution algorithm, and is insensitive to the ma- 
trix element accuracy (i.e. ANk). 

b. Energy The field energy is evaluated according to en- 
ergy functional Eq. (13). Unlike normalization, which can be 
calculated to numerical precision, the energy is Umited to the 
precision with which we can evaluate the dipole energy. For 
the results in Table III the energy functional is evaluated for 
the same value of ANj. as was used for the evolution under 
consideration. These results reveal a similar convergence be- 
havior to that observed for 5N. 



c. Angular momentum For the dipolar system the 
anisotropic nature of the long-range interaction leads to inter- 
esting dynamics of the angular momentum, which we discuss 
further in Sec. VG2. However, for the case of a spherical 
trap the z component of angular momentum is conserved. To 
characterize this we evaluate 



j d3a;V'^(x,i)L,Vc(x,i), (76) 



where is the z component of L = — zfi.x x V. Like normal- 
ization (and in contrast to the energy), the angular momentum 
can be evaluated efficiently and to numerical precision using 
the step operator formalism, as discussed in Ref. [40]. The 
results in Table III show that SL^ appears to converge quite 
slowly in AA^^, and is conserved at the 10~^ level for our 
ANk = 40 simulations. This may indicate an important con- 
sideration for the dipolar PGPE, and we discuss this further 
below. 

d. Field convergence The quantity SX indicates the ex- 
tent to which the field evolution has converged. The results 
for SX in Table HI have been computed by comparing each 
case to a more accurate calculation with a relative tolerance 
of 10 ~^ and the same AN/, value. These results are insen- 
sitive to ANk and show rapid convergence as the evolution 
tolerance is decreased. However, an important dependence on 
ANk is revealed by computing SX', defined as in Eq. (75), but 
by comparing the against a c'^(T) for a different (i.e. larger) 
ANk value. These results, presented in Table III for the case 
where the accurate solution uses a relative tolerance of 10~® 
and ANk = 50, reveal a much slower convergence in the 
parameter AA^^, with a very weak dependence on evolution 
tolerance. This appears to be due to the rather slow conver- 
gence of the high energy matrix elements with ANk as noted 
earlier (e.g. see Sec. VD). These results serve to illustrate an 
important point: Our algorithm in a fixed ANk subspace is 
well-defined and displays good convergence primarily depen- 
dent on the evolution tolerance. 

We note that the individual simulations reported in Table III 
took between 3 minutes (~ 350 steps with ANk = 0) and 2 
hours (~ 2000 steps with AA^^ = 40) using unoptinuzed sin- 
gle CPU code running on a shared cluster of 2.66GHz Clover- 
town Xeons. 



G. Convergence of thermodynamic predictions 

The important question we have yet to address is: what ac- 
curacy is required to perform a useful PGPE simulation? In 
general the answer to this question will depend on the partic- 
ular application of interest, and in this final part of the paper 
we will present some illustrative examples. 

For deterministic applications, such as solving a T = 
Gross-Pitaevskii equation from a well-defined initial state, the 
small errors in the matrix elements will cause errors to ac- 
cumulate leading to a practical time limit for the duration 
over which a calculation can be considered to be reliable. In 
contrast, the PGPE theory is typically operated in an ergodic 
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regime of evolution, in which we only aim to specify or mea- 
sure macroscopic features of the field. An approximate treat- 
ment of the dipole interaction (e.g. all matrix elements at the 
10 ~^ level of accuracy or better) would seem to be more than 
adequate for such applications, as long as our approach does 
not break important symmetries of the system, e.g. allowing 
constants of motion to change appreciably with time so that 
the system relaxes to the wrong equilibrium state. 

To investigate these issues we simulate the evolution of the 
dipolar PGPE in a finite temperature regime, and explore how 
changing AN), affects its predictions. To do this we prepare 
a random state of energy E = 10.0, for an isotropic harmonic 
trap with C = 500, D = 500 and 

Ccut — 23. We use this 
state as the initial condition for 8 simulations which differ in 
AN), from to 28. In each case we propagate the dipolar 
PGPE, using the adaptive step Runge-Kutta algorithm with a 
tolerance of 10"'', for T = SOtt (i.e. 40 trap periods), saving 
the field at 1600 equally spaced times during the evolution. 



1. System width 

The randomly generated initial state used in the PGPE is 
an atypical (far from equilibrium state) and will evolve for 
some initial period until the system explores more typical 
microstates (i.e. rethermalizes). After this initial period we 
can compute ensemble averages of equilibrium parameters by 
making use of the system's ergodicity. 

A simple macroscopic parameter to compute is the mean 
system width, as characterized by the position variance, e.g. 
Wx{t) = {x'^{t)) - {x{t))'^ in the x direction, where 



(a;"(t))= / d'xx"\M^,t)\' 



(77) 



is the instantaneous moment. To make equilibrium predictions 
it is useful to calculate the averaged width, which we calculate 
using time-averaging, i.e. the time averaged moment is given 
by 



_ 1 ^= 



(78) 



where Ns is the number of samples used. We avoid writing 
the similar expressions for Wy and Wz- 

In what follows we let the system thermalize for the first 
10 trap periods (in practice most large scale motion damps in 
the first few trap periods), and then perform time averaging 
using Ns ~ 1300 states over the subsequent 32 trap period 
evolution. 

The results for the position width are shown in Fig. 5. In- 
terestingly the width of the system in the z direction is greater 
than the x and y directions even though the system is in an 
isotropic harmonic trap. This asymmetry arises from the po- 
larization of the dipoles in the z direction which causes the 
system to slightly elongate to reduce the dipolar interaction 
energy. We note that there is no clear change in the results 
with ANk [53] and the improved accuracy associated with in- 
creasing ANk is clearly unimportant in this case. The states 
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FIG. 5: (a) Evolution of z width for simulation with ANk = 
(grey line) and AA^fe = 28 (black line), (b) Dependence of the time 
averaged predictions for the variances of a dipolar Bose gas on ANk. 
Results for Wx (circles), Wy (triangles), and Wz (squares) in units 
of the harmonic oscillator length squared (xq = h/nuo). Shaded 
region in (a) indicates the states used for time averaging. The shaded 
regions in (b) characterize the spread in results. AU results computed 
using the corrected dipole interaction. 



over which time averaging is performed are indicated by the 
shaded region in Fig. 5(a). Interestingly the breadth of this 
region (chosen to match the range of the equilibrium width 
dynamics) is 20 times larger than the shaded region shown in 
Fig. 5(b) to indicate the spread in the averaged z width re- 
sults. This suggests that while the width dynamics are quite 
appreciable, the averages are very well-defined. Longer time 
averages could be used to further refine these predictions. We 
also note that the larger variation in the x and y variances seem 
to result from strong collective dynamics associated with the 
non-conservation of angular momentum, which we discuss 
below. 



2. Angular momentum evolution 

The anisotropic (non-central) nature of the dipole interac- 
tion means that angular momentum is not conserved even for 
the case of a spherical external potential. Indeed, as can be 
shown (see Appendix B) the evolution of the angular momen- 
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turn is given by 

dt 

d{Ly) 

dt 
d{L,) 

dt 



-A-kD j <fkh{\<i)h{-\<i) 
4ttD j d^kh{k)h{-k) 



kykz 



kz kx 



0, 



(79) 



(80) 



(81) 



(for the isotropic trap case), revealing that the invariance of 
rotations about the polarization direction leads to conservation 
of the z component of L. This motivated the definition of SL^ 
as a numerical check in Sec. V F. 



-0.5 



-0.05 




0.05 



¥IG. 6: Angular momentum evolution, (a) x and y components of 
angular momentum in units of h over a small time segment of the 
simulation for AiY^ = 0. (b) The z component of angular momen- 
tum of the whole evolution for simulations of various /S.Nk values. 
Parameters are the same as for the results presented in Fig. 5 and 
time is measured in units of inverse trap frequency (l/w). All results 
computed using the corrected dipole interaction. 

In Fig. 6(a) we show the evolution of the x and y compo- 
nents of angular momentum for our dipolar simulations. As 
suggested by Eqs. (79)-(81), the x and y components of an- 
gular momentum show strong dynamics. These dynamics are 
a contributing factor to the slightly larger spread in results for 
the position variance in the x and y directions relative to the z 
direction, as seen in Fig. 5(b). 

In Fig. 6(b) we examine the evolution of L^. According to 
Eq. (81) Lz should be conserved, and so the dynamics of this 
quantity indicates inaccuracy in our algorithm. The various 



curves in Fig. 6(b) indicate that as ATV^ increases, the drift 
in Lz decreases. In some applications of the dipolar PGPE 
theory, e.g. in studies of vortices, careful attention to con- 
servation will be prudent and will demand the use of a large 
AA^fc. However, for many applications the quasi- stationary 
behavior of observed in the AiVfc = case will be ade- 
quate to make reliable predictions (e.g. our position variance 
results appear insensitive to AN).). 



VI. CONCLUSIONS 

In this paper we have presented a numerical method that 
allows us to extend the PGPE theory to include long-range 
dipolar interactions. We have used a range of tests to charac- 
terize the numerical accuracy of our scheme and the conver- 
gence with increasing order of fc-space quadrature grid. These 
results show that use of the corrected dipole potential is a sig- 
nificant improvement, and that our approach is sufficiently ac- 
curate to make reliable physical predictions in the context of 
finite temperature c-field calculations. Many aspects of the 
formalism we have developed are quite general and would 
easily allow us to apply the method to a wider class of long- 
range interactions. 
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APPENDIX A: RANDOMIZED STATE 

We generate a pseudo random slate based on a linear con- 
gmential generator, with recurrence relation 



Xn+\ = mod m{aXn + c), 



(Al) 



with a = 16807, c = 0, and m = 2^^ - 1. 

We prepare a set of complex random numbers defining the 
classical field, Cap-y. To do this we map the quantum number 
tuples {a, /?, 7} to a unique integer value, n, according to 



n = a + M^P + MJ7. 
We then specify our classical field state as 



(A2) 



(A3) 



where Xn' is the sequence generated by (Al) with seed 
Xq^^ = 10^, and Xn ^ is the sequence generated by (Al) with 
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= 10^ Thus, we have 



Co = cooo = 0.0466 + 0.4657i, 
ci = cioo = 0.6369 + 0.3693i, 
C2 = C200 = 0.8143 + 0.1432i, 



(A4) 
(A5) 
(A6) 
(A7) 



APPENDIX B: ANGULAR MOMENTUM EHRENFEST 
RELATION 

Given angular momentum operator L = —itix. x V, the 
standard Ehrenfest result for the GPE angular momentum is 
given by 



ih^ = {L[Vo + n. 



(Bl) 



The effect of the harmonic trap potential, Vq, on the angular 
momentum evolution is well-understood. Here we will fo- 
cus on the case of an isotropic trap (i.e. [L, Vq] = 0) so that 
the evolution arises from the effective dipole potential, i.e. 
ihd{L)/dt = (L^). 
Taking the Fourier transformed form of 5> (see Eq. (12)) 

(L$) = j d^xijj^iix) j d3fcVD(k)n(k) (Le*''''') Vc(x). 

(B2) 

We can use the self-duaUty of angular momentum operators 
under Fourier transform, that is 



j d^k-Le^^ '' = - J d^k 



Le 



ik-x 



(B3) 



where L is the representation of angular momentum in k- 
space, i.e. = —ih{kxdk — kydk^). We then find 



(L$) 



-I 



so that 



ih 



j d^kh{-\<i)t (yD(k)n(k)) , 

d^A;n(-k)n(k)LVD(k) 
+ j d3fc^yzj(k)L(ri(-k)n(k)), 
d3fcn(-k)n(k)LVD(k), 



,3, n(k)n(-k)~--_ 



= / 



di , d^^-'^^^^Vniy)- 



(B4) 
(B5) 

(B6) 
(B7) 



We now make use of the Cartesian components of L in spher- 
ical co-ordinates: 



coscpk d . d 
tan Ok 0(pk otik 
sm(j)k d 6 ^ 

tan Ok dcpk d9k ' 



= —ih- 



d 



(B8) 
(B9) 
(BIO) 



where (pk is the azimuthal angle from kx in the kx-ky plane. 
For the dipolar potential, we find 



" ^ ' |k|' ■ 

k 



and the angular momentimi equations 



d{Lx) 
dt 

d{Ly) 
dt 

dt 



-47r£» j d^kn{k)n{-k) 

AttD j d''kh{k)fi{-k)^ 
0, 



kykz 

2"' 



(BID 

(B12) 
(B13) 

(B14) 
(B15) 
(B16) 



which should provide useful consistency conditions for nu- 
merical simulations. d{Lz)/dt = is expected from 
the cyhndrical symmetry of Fd(x) about the polarization 
axis. We can also see from Eq. (B14) that if n(k) = 
h{kx,\ky\,\k^_\) i.e. is reflection symmetric in the ky and 
kz directions, then d{Lx)/dlt = 0. Similarly, if n(k) = 
n{\kx\,ky, \kz\), then d{Ly)/dt = 0. We note that n(— k) = 
n(k) holds when n(— x) = n(x) so that eigenstates of parity 
will conserve L. Consequently the evolution of a spherically 
symmetric state into a cylindrically symmetric state should 
conserve angular momentum. We have not included boundary 
terms in this derivation which arise from the projector, and fu- 
ture work will be to assess at what level they may contribute 
(e.g. see [54]) 



APPENDIX C: ANALYTIC EVALUATION OF THE PURE 
DIPOLE MATRIX ELEMENTS 

In this appendix we derive an analytical expression for the 
pure dipolar matrix elements, as given by Eq. (66): 

Zi% = Ci%jd\d\'e-(^'^y'^^'^Hi%{x) 

X y^(x-x')e-(-"+''"+^")if:^:(x'), (CI) 



where 



and 



(C2) 



H'^^^ix) = Hs{x)H,{v)Hc{z)H^{x)Hp{y)H^{z). (C3) 
Using the convolution theorem with 

AttD ( 3fc2 N 



^{Vd(x-x')} = 



-1 , (C4) 
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and 



(C5) 

where 

the dipolar matrix elements can be evaluated from 

3A^^ 



a/57 ^ c«P7 I ^ ^ 



2 a/37 



1 2 I 1,2 I 1,2 



i/3 



\a-5\ _ 
Min[a,i5] I 2 ^ "^Min[/3, 



-1/3-^1 



-I7-CI 

^Min[7,C] 



where 

RiSeC _ riSeC ^ / -i \Max[a,i5]+Max[/3,e]+Max[7,C] 

-°a/37 ~ *-^a/37'-^a^7^ -^-f 

^ ^a+/3+7+5+e+C2-f+Min[a,5]+Min[/3,e]+Min[7,f] 

X Min[a,5]!Min[/3,e]!Min[7,C]!. (C8) 



Expressing the associated Laguerre polynomials in terms of a 
finite sum: 



1) LUx) = J2i-iy 



{n + ky. 



m=0 



(n — m)!(fc + m)!m! 



(C9) 



we find that if a — (5 or /? — e or 7 — ^ is odd then Z^'^j^ = 0. 
(CV^hen a — S and (3 — e and 7 — C are even Eq. (C7) reduces to 



''a/37 



a /? 7 Min[Q.5] Min[/3,€] Min[7.C] 

E E E 

jl =0^2=0^3=0 j4=0 j5=0 j6=0 



2j25 + l/3-e|-2 

n (2(ji2+j45) + |a-<5| + |/3-e|-j7) 

J7=0,2,4... 



(_l)i2i(i-2j-i4+l/3-^l+l7-CI) (2j _ 6i36 + \a-6\ + \p-e\- 2|7 - C|) (-1 + 2ju + \a - 5|)!! (-1 + 2j25 + |/3 - e|)!! 



(ii!j2!i3!) J4!i5!i6!(a - jiV-iP - i2)!(7 - J3)!04 + \a- (5|)!(i5 + 1/3 - e|)!(i6 + |7 - CI)! 

r [1 + j - he +\a- S\/2 + 1/3 - e\/2] T [j36 + (1 + I7 - CI)/2] 

(ji4 + |a - <5|/2)! (3 + 2j + |a - (51 + 1/3 - e| + |7 - CI) (Min[a, S] - j4)!(Min[/3, e] - j5)!(Min[7, C] - je)! 

I 



(CIO) 



where 



Bi%a\f3\^\{Mm[a,5] + \a-6\)\ 
X (Min[/3, e] + 1/3 - e|)!(Min[7, C] + I7 - CI)!, 

(Cll) 



j = h + h + h + k + h + k and jab = 3 a + k- 



Equation (CIO) can be readily evaluated and serves a direct 
comparison for the numerical integration. 
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